Instanton filtering for the stochastic Burgers equation 
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We address the question whether one can identify instantons in direct numerical simulations of 
the stochastically driven Burgers equation. For this purpose, we first solve the instanton equations 
using the Chernykh-Stepanov method [Phys. Rev. E 64, 026306 (2001)]. These results are then 
compared to direct numerical simulations by introducing a filtering technique to extract prescribed 
rare events from massive data sets of realizations. Using this approach we can extract the entire 
time history of the instanton evolution which allows us to identify the different phases predicted by 
the direct method of Chernykh and Stepanov with remarkable agreement. 
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Introduction Understanding intermittency in turbulent 
flows is still one of the open problems in classical physics. 
More than 15 years ago, for certain systems like the prob- 
lem of passive advection and Burgers turbulence the door 
for attacking this issue was opened by getting access to 
the probability density function to rare and strong fluctu- 
ations by the instanton approach [ll-Q ■ In this letter we 
concentrate on rare fluctuations in Burgers turbulence. 
For that case, Gurarie and Migdal Q introduced the in- 
stanton approach and were able to calculate the instanton 
contribution to the right tail of the velocity increment 
probability distribution function (PDF). In succeeding 
work, Balkovsky et al. Q were able to characterize the 
left tail of the increment PDF making use of the Cole- 
Hopf transformation 0, Q . These analytical results were 
confirmed by direct numerical solution of the instanton 
equations by Chernykh and Stepanov Q ■ 

The open question remained whether one can observe 
or identify the instanton in numerical simulations of the 
stochastic Burgers equation. The answer is not obvious, 
since one could argue that, perhaps, the contribution of 
the instanton is exponentially small such that instantons 
are only relevant to such rare events that they are not in- 
teresting from a practical point of view at finite Reynolds 
numbers. In this Letter, however, we find that already at 
moderate Reynolds numbers the instanton can be identi- 
fied in data sets of simulations of the stochastic Burgers 
equation. This gives a positive answer to this important 
question. In particular, we show by introducing a partic- 
ular filtering technique that all phases of the instanton 
evolution can be recovered from data sets of simulations 
of the stochastic Burgers equation. 

The Letter is organized as follows: we first review 
the path integral formulation for Burgers turbulence 
and revisit the algorithm introduced by Chernykh and 
Stepanov to solve directly the instanton equations. 
We then describe our numerical simulations to obtain 
sufficient statistical data (« 10 7 realizations of complete 
stochastic Burgers simulations using CUDA graphics 



cards) necessary for our instanton filtering. Using this 
enormous data set we apply our instanton filtering pro- 
cedure and compare the results with the direct instanton 
simulations. A conclusion and outlook summarize the 
Letter. 

Action functional and instanton equations. We consider 
the stochastically driven Burgers equation given by 
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with a noise field (j) that is 5-correlated in time and has 
finite correlation in space with correlation length L, more 
precisely 



{<t>(x,t)cl>{x',t')) 

xO) 



6(t-t') X ((x-x')/L), (2) 
(1 - z 2 )e- x2/2 ■ (3) 



While the precise form of x is n °t important for the re- 
sults of our work, we chose this particular form in order 
to have the same setup of the problem as in previous 
studies by Chernykh and Stepanov. Using the functional 
path integral introduced by Martin-Siggia-Rosc/Jansscn- 
dc Dominicis [§-fll|. the PDF of the velocity gradients 
u x (t = 0, x = 0) is written as 



V(a) = JvuVpdFexp (-S(u,p, F) 



(4) 



with the action S given by 

S = ^ J dt J dxi dx 2 p(xi,t)x(xi - x 2 )p(x 2 ,t) 
-i dt dx p(u t + uu x - vu xx ) 

J — oo J 

-Av 2 Ti{u x {0,Q) - a) , (5) 

where T results from the Fourier transform of the 5- 
function for the observable u x (t = 0, x = 0) = a. The 
saddle point (instanton) equations for the fields (u,p) 
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yielding the largest contribution to the path integral for 
strong gradients are then given by 



hi 



p t + up x + vp x 



x')p(x',t)dx' (6a) 
Aiv 2 F5{t)5'{x) . (6b) 



The Chernykh-Stepanov algorithm revisited. The algo- 
rithm proposed by Chernykh and Stepanov for solving 
the system of partial differential equations for the fields 
u and p can be summarized as follows: the diffusion 
terms in the equations ([6]) define the temporal direction 
of the numerical integration of the equations, meaning 
that u is integrated forward in time while p is integrated 
backwards. The right-hand side of equation ()6b[) poses 
the initial condition p(t = Q,x) = —Aii/ 2 J-6'(x) and the 
starting step is obtained by setting u(t, x) = 0. Equation 
(|6b|) is then solved backward in time up to a large 
negative time mimicking -co. The obtained solution 
p(t,x) is used in the right-hand side of equation (|6a|) 
such that this equation can be solved forward in time. 
This procedure is then iterated until convergence. For 
higher gradients, a stabilization of this iteration has to 
be applied, details of which can be found in 0. While 
Chernykh and Stepanov use a stabilized finite difference 
scheme with an implicit first-order time integration, 
we utilize a second-order Adams-Bashforth temporal 
integration for a pseudo-spectral method. We also note 
the similarity of the system ^ to equations that arise in 
the context of transition probabilities [12, 13 1. Although 
the boundary conditions are different, the above system 
of instanton equations can, in principle, also be solved 
numerically by minimizing the corresponding action 
using a L-BFGS scheme. We found, however, for 
the case under consideration, the propagation-based 
Chernykh-Stepanov scheme numerically much more 
efficient. Therefore, the Chernykh-Stepanov scheme 
might be an interesting alternative to compute transition 
probabilities. A detailed comparison of both schemes 
is beyond the particular scope of this paper and will 
presented elsewhere. 

Parallel simulation of the stochastic Burgers equation. In 
order to generate data from simulations of the stochas- 
tic Burgers equation, we need to solve eq. (U} with the 
appropriate right-hand side. For the generation of the 
stochastic force field <f>, at each step in time, we draw a 
vector r of appropriately scaled normally distributed ran- 
dom numbers. The size of the vector corresponds to the 
discretization in x. This vector r is then multiplied by a 
matrix A resulting from the Cholesky decomposition of 
the (discretized) correlation matrix C. Note that naive 
discretization of x may lead to a C that, due to finite 
machine-precision, is not positive-semidefinite. There- 
fore we used the algorithm introduced in [14J in order to 
obtain a matrix C that is positive-semidefinite and suffi- 
ciently close to C. Note that this method of generating 
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FIG. 1: Comparison of the filtered velocity field (u s hift e d(£, x)) 
(top) and the instanton field u(t, x) (bottom) at times t = 
(solid) and t = —1.75 (dashed). 



the noise term in eq. (JXJ) is different from Fourier-based 
methods that are commonly used for such simulations 
15l[l6| and more computationally expensive. The reason 
for the presented choice was motivated by the necessity 
to numerically generate noise that closely imitates the 
fluctuations assumed in the instanton analysis. 

As the size of an individual realization is small enough 
to fit on a single graphics card with its complete 
history, the whole simulation is performed in CUDA 
alone. Matrix- Vector-operations are realized using the 
cuBLAS-packagc, the fast Fourier-transform is provided 
by cuFFT. Since the filtering and shifting procedure was 
also implemented in CUDA, a whole bulk of simulations 
is performed and filtered independently on the GPU. Av- 
eraging over different CUDA-processes occurs after com- 
pletion of a bulk and is performed via MPI. Thus, both 
expensive device-to-host copies and high-latency network 
communication are minimized. Because of the stochastic 
independence of realizations, this method scales linearly 
with the number of graphics cards. 

Extracting the instanton. In order to provide a sufficient 
data set for the extraction of the instanton from simu- 
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Run 1 


1024 


0.039 


0.406 


1 


40 


0.3 


4.586 


0.99 


10.5 


Run 2 


1024 


0.039 


0.464 


1 


40 


0.38 


2.691 


0.97 


0.410 


Run 3 


1024 


0.039 


0.481 


1 


40 


0.41 


2.33 


0.95 


0.052 



TABLE I: Parameters of the numerical simulations. JV: num- 
ber of collocation points, dx: grid-spacing, 77 = (u 3 /e^) 1 ^: 
Kolmogorov dissipation length scale, L: correlation length 
of forcing, Lbo X : domain length, v: kinematic viscosity, £k: 
mean kinetic energy dissipation rate, Tz, = L/u rm8 : large- 
eddy turnover time, #hits (%): percentage of hits with pre- 
scribed velocity derivative. 
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FIG. 2: Comparison of the filtered velocity field (u s hifted(i, a;)) 
(top) and the instanton field u(t, x) (bottom) as a space-time 
contour plot. 



lations of the stochastic Burgers equation, wc conducted 
the following numerical experiment: We started the in- 
tegration of the stochastic Burgers equation from zero 
initial conditions for the velocity field u(t = t m [ n ,x) = 
at a large negative time i m j n up to the final time t = 0. 
Typical parameters are summarized in table [TJ The ini- 
tial time £ m i n was chosen corresponding to the instanton 
equations (|6]) and consists of more than 10 integral times 
Ti,. This single experiment was repeated « 10 7 times us- 
ing the 64 CUDA Tcsla 1060 graphics on the Bochum 



GPU Cluster and the 96 CUDA Fermi 2050 graphics 
cards on the CUNY GPU Cluster. The total simula- 
tion length obtained by this parallelism corresponds to 
« 10 s integral times Tl- Each of these simulations was 
analyzed in the following way: We prescribed a small 
interval around a given value of the velocity gradient 
u x {t = 0,x) = a at the final time t = and searched 
for the maximum velocity gradient in the numerical so- 
lution at that time. If we find that the maximum veloc- 
ity gradient falls into the desired interval, we shift the 
field in space such that the location of the maximum 
velocity gradient is located at x = 0. In addition, we 
also shift the forcing field <f>(t, x) in the same way. The 
averaging procedure now consists of taking the average 
of all those shifted fields w s hiftcd(i 7 x) and Shifted {t, x). 
We thus obtain an ensemble average (u s hifted(ij x )) an d 
(0shifted(£> x)) in space and time. Since the forcing field is 
(5-correlated in time, it is obvious that in order to extract 
information of the averaged forcing field an enormous 
number of realizations is necessary. This numerical pro- 
cedure now complies with the path integral formulation 
for the observable 0(u) — (6(u x (0, 0) = a)) 

(0(u)) =[vf 0{u)8{u t + uu x - vu xx - 0)e-<*'X -1 *)/ 2 

(7) 

which is the starting point for the Martin- Siggia-Rose 
formulation. Thus, for sufficiently strong velocity gra- 
dients u x (0,0) = a, the important question and con- 
jecture is whether the averaged solutions (ushifted(i) %)) 
and {(/'shifted (t, x)} coincide with the instanton solution 
of ([6]). Especially if this conjecture is true, then the av- 
eraged optimal force (4>shifted(t,x)) should coincide with 
the right-hand side of equation (|6a|) 

(Shifted = -i J x(x - x')p(x' ,t)dx' (8) 

where the auxiliary field p(t, x) is obtained from the di- 
rect Chernykh-Stepanov algorithm. Fig. ([1]) shows the 
filtered field (u s hiftcd(t, x)) (top) and the instanton field 
u{t,x) (bottom) at the final time t = and at an earlier 
time t = —1.75 showing the instanton in a different phase 
(see also the sketch of the instanton phases in Fig. 8 in 
[7|). The agreement is remarkable. Especially the center 
region is precisely reproduced by the stochastic simula- 
tion, while the sides are less pronounced. In order to get 
a complete overview of the time history of the instanton 
and the filtered field, Fig. ([2} depicts a contour plot of 
the whole space-time domain. Although the filtered field 
shows a slightly shorter extent in time, the congruence is 
clearly visible. 

The rareness of the filtered events has a strong impact 
on the agreement between the instanton approximation 
and the full stochastic simulation. In order to demon- 
strate the varying resemblance to the instanton approxi- 
mation, we alter the probability of reaching a prescribed 
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FIG. 3: Comparison of the instanton field u(t, x) (solid) to 
stochastic simulations with varying hit percentages (sa 10% 
(dashes), w 0.5% (dots), « 0.05% (small dots)) for t = 
(top) and t — —1.75 (bottom). Agreement with the instanton 
approach increases with decreasing hit percentage. 



velocity gradient by changing the kinematic viscosity v. 
Fig. (0} shows the filtered field (u s hihed(t, x)} and the in- 
stanton field u(t,x) for three different hit percentages. 
As the rareness of the event increases, accordance with 
the instanton grows considerably. Especially the velocity 
gradient in the origin is only reproduced when the events 
are rare. Notably this effect docs not depend on the 
Reynolds number or shock strength, but on the scarcity 
of the event alone. 

An additional feature of this filtering approach is 
that not only the instanton velocity field could be 
extracted but also the time history of the auxiliary 
field p(t, x) and of the optimal force field. At time 
t = 0, the auxiliary field is given by its initial condition 
p(t = Q,x) = —Aiv 2 F5'(x) and produces the force 
term Av 2 !F'x'(x) on the right-hand side of eqn. (|6a[). 
A comparison of this term with the filtered force field 
(<Ashifted(i, x)) is depicted in Fig. (j4]), which shows a 
remarkable agreement. 
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FIG. 4: The filtered force field ((^shifted (t, x)} (dashed) and 
the analytical force field Av 2 \Fx' ( x ) (solid) at time t = 0. 



Conclusions and Outlook In this Letter we studied the 
question whether the instanton solution for Burgers 
turbulence is "real" , e.g. can be observed in stochas- 
tically driven simulations. The positive answer to this 
question is remarkable since we observe the instanton 
for a moderate Reynolds number and thus for moderate 
(and not extreme) values of the velocity gradient. In 
principle, this filtering method also allows a further 
study of moderately scarce events to determine where 
fluctuations around the instanton appear, how they look 
like and how they modify the action integral and thus the 
PDF. Our findings could also open the door to the issue 
why the PDF for the very left tail of velocity gradients 
could not be observed in high resolution numerical exper- 
iments of Gotoh [3] . Although this is out of the scope 
of the present Letter, work in this direction is in progress. 
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